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OO Abstract. Betelgeuse, the bright red supergiant (RSG) in Orion, is a 

runaway star. Its supersonic motion through the interstellar medium 
_ has resulted in the formation of a bow shock, a cometary structure 

pointing in the direction of motion. We present the first 3D hydrody- 
namic simulations of the formation and evolution of Betelgeuse's bow 
shock. We show that the bow shock morphology depends substantially 
on the growth timescale for Ray leigh- Taylor versus Kelvin-Helmholtz 
Oh instabilities. We discuss our models in light of the recent Herschel, 

^ GALEX and VLA observations. If the mass in the bow shock shell 

^ is low ('^fewx 10~^ M©), as seems to be implied by the AKARI and 

Herschel observations, then Betelgeuse's bow shock is very young and 
is unlikely to have reached a steady state. The circular, smooth bow 
shock shell is consistent with this conclusion. We further discuss the 
I implications of our results, in particular, the possibility that Betelgeuse 

^ may have only recently entered the RSG phase. 

o 
o 

04 1 Introduction 

^T) From the shoulder of Orion 'The Hunter' (Greek mythology), to stories about 
the fierce, red lion (Southern African mythology), Betelgeuse has long been a 
prominent part of the night sky. As the nearest and brightest star of its kind, 
\ I Betelgeuse is now considered the prototype red supergiant (RSG). Estimates of its 
^ mass range from 8Mq - 20 and although it is very cool (Teff~3 300 K), it is 
highly luminous (L* ^lO^L©) due to its large stellar radius {R^ ^lOOOR©) (see 
for example. Smith et al. I2009| Neilson et al. I201ip . Its tenuous atmosphere is only 
loosely bound, consequently it loses ~ 2 — 4 x 10~^ M0yr~^ via a slow, ~17kms~^ 
wind (Noriega-Grespo et al. I1997[ Bernat et al. .1979 ). The mechanism by which 
this material is lost is still unclear, but the process has occurred for thousands of 
years forming an extensive circumstellar envelope (GSE) of gas and dust. 
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Noriega- Crespo et al (|1997p successfully imaged the CSE at 60 and 100 pm 
using IRAS. They detected a bow shock arc ~6 arcmin in radius to the north-east 
of Betelgeuse and a mysterious linear 'bar-like' structure, located just ahead of the 
arc. A decade later, Ueta et al. (|2008|) confirmed the detection, imaging part of 
the bow shock arc and bar in the far-infrared with AKARI. More recently, high 
resolution Herschel observations revealed that the bow shock consists of multiple 
arcs (Cox et al. [2QT21 Decin et al. I2QT2]) . Le Bertre et al. (f2QT2)) identified a faint, 
0^4 LE'X far- ultraviolet arc at the same position as the outermost Hershel one. The 
bar, however, was not detected at these shorter wavelengths. With the VLA^ they 
also found atomic hydrogen coincident with the bow shock and an inner, comet ary 
shaped detached shell of HI emission ~4 arcmin in diameter. 

Betelgeuse is moving supersonically relative to the local interstellar medium 
(ISMQ and its bow shock is formed by the collision of its stellar wind with this 
medium. Assuming it has reached a steady state, the bow shock can be used to 
probe the physical properties of these interacting flows. The bow shock 'radius', 
known as the stand-off distance, i^so, is the location where the ram pressures of 
the ISM and stellar wind are in equilibrium, and is given by: 

PiSM^* = MwVw/47ri?so ' (1-1) 

(assuming a spherical wind) where is the wind mass-loss rate, pisM and are 
the density of the ISM and stellar wind, respectively; is the velocity of the star 
with respect to the ISM, and is the stellar wind velocity. Assuming momentum 
conservation and that the stellar wind and ISM mix and cool instantaneously (the 
thin-shell approximation), the shape of the bow shock is given by: 

R{0) = i?so cosec 6> ^3(1 -6> cot i9) , (1.2) 

where 9 is the polar angle measured from the axis of symmetry (Wilkin [19^961) . 

Utilising these analytic models and current estimates for Betelgeuse's wind and 
distance, Ueta et al. (2008) derived a space velocity of = 40 Hj^^^^ kms~^ with 
respect to the local ISM. Estimates of the ISM density, uh, range from 0.3 cm~^ to 
1.5 - 1.9 cm~^, thus Betelgeuse's space velocity is likely to be between 73 kms~^ 
and 28 kms~^, respectively. Mohamed et al. (2012) simulated models for these 
parameters and compared the results to the IRAS and AKARI observations. In 
this paper, we highlight the main points of that study and discuss the conclusions 
in light of the recent Herschel^ GALEX and VLA observations. 

2 Model 

The bow shock is modeled in 3D using Smoothed Particle Hydrodynamics (SPH), 
a Lagrangian method particularly suited to studying hydrodynamical flows with 



-•^ Although Betelgeuse is the only RSG with a bow shock, theoretical models predict that up 
to 30% of RSGs can be runaway stars (Eldrige et al. I2011|) . The origin of the Betelgeuse's 
high space velocity is unclear but may be due to a dynamical ejection from a cluster and/or a 
supernova kick. 
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Fig. 1. Evolution of the ratio of i?(0°)/i?(90°) for a slow, = 32kms~^ model (ma- 
genta), and fast, — 73kms~^ model (dashed blue), compared to the analytic, steady 
state value, (solid black line). The observed Herschel/ AKARI ratio is ^0.7 (ar- 

rows) . 



arbitrary geometries. Throughout this study we use the GADGET-2 SPH code 
( [Springel 2005) which we have modified to include stellar winds (Mohamed & 
Podsiadlowski [2007 ) , an ISM flow (Mohamed I2010p , and atomic and molecular 
radiative cooling (Smith & Rosen l2003p . 

For numerical convenience, we select the stellar rest frame with the star located 
at the origin of a rectangular box {x^y^z = 0, 0, 0). Given the uncertainty in 
Betelgeuse's mass-loss mechanism, we do not model the wind acceleration in detail. 
Instead wind particles are injected isotropically at a radius, i^inner'^lO^^ cm, with 
velocity, 'U^^lTkms"^, and temperature T^^IOOO K. The result is a smooth, 
constant outflow of material at a rate of 3.1 x lO~^M0yr~^. The ISM is also 
assumed to be homogeneous and flows past the star in the direction of the x 
axis, interacting with the stellar wind as it does so. We model a range of ISM 
densities, uh = 0.3, 1.0, 1.5, and 1.9 cm~^, with corresponding stellar velocities, 
V:^ = 73, 40, 32, and 29 kms~^, respectively. These number densities lie at the 
boundary between typical values expected for either a warm or cold neutral ISM, 
so we assume temperatures based on the phase diagram of the standard model of 
Wolfire et al fT995), e.g., their Fig. 3d. The temperatures, Tjsm, are 8 000, 1600, 
1 000, and 650 K, respectively. Additional models are also run to investigate the 
effect of varying the ISM temperature, degree of cooling and numerical resolution. 

The numerical method and model set up were tested with an adiabatic model. 
The results were consistent with both theoretical expectations and previous studies 
(e.g., Wilkin [19961 Brighenti & D'Er cole 1995^. 
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Fig. 

kms" 



2. Hydrogen column density for the slowest (29 kms~ ) [a] and the fastest (73 
[b] models after 20 000 years (left) and 32 000 years (right), respectively. For 
animations of the bow shock evolution, see Mohamed et al. l20T2l 



3 Results 



The simulations begin at the start of the RSG phase. As the stellar wind collides 
with the ISM, material accumulates at the contact discontinuity, where part of 
the kinetic energy of the gas is thermalised. The heated ISM and stellar wind 
expand outwards from either side of the contact discontinuity; the former pushes 
into the ISM, this is known as the forward shock, and the latter pushes into the 
stellar wind, known as the reverse shock. Although the stellar outflow is initially 
spherical, it becomes increasingly parabolic as the star moves through the ISM. 
Eventually a steady state is achieved at which point the global morphology is 
described by Eq. |1.2| From the models we derive the ratio of the bow shock radius 
at angles l9 = 0° and l9 = 90°, i?(0°)/i?(90°), as a function of time (shown in 
Fig. [T]). It takes several dynamical timescales for the bow shock to achieve the 
equilibrium value i?(0°)/i?(90°)=l/V^. 

Although all the models exhibit a similar global structure, the flow character- 
istics on smaller scales differ considerably due to the growth of Ray leigh- Taylor 
(R-T) and Kelvin- Helmholtz (K-H) instabilitie^ (see Fig. [2|. In the 'slow' mod- 
els, nn ^lcm~^ and <40kms~^, the strong cooling reduces the thermal pres- 
sure of the gas enabling further compression in the forward and reverse shocks. 
The greater post-shock densities reduce the growth timescale for R-T instabilities. 
This, along with the slow space motion producing less shear, causes the R-T 'fin- 
gers' to develop faster than the K-H 'rolls'. These bow shocks consist of a thin, 
smooth outer shock and a contact discontinuity that is highly distorted by R-T 
'fingers'. In the column density plots, the small-scale R-T instabilities result in a 
clumpy, knot-like sub-structure that becomes one of the dominant features of the 
bow shock, particularly when viewed at large inclination angles (see Fig. [s] [top]). 



^ Ray leigh- Taylor instabilities are 'finger-like' protrusions that occur when a light fluid is 
accelerated into a denser fluid. Kelvin-Helmholtz 'rolls' or 'eyes' are excited by the shear produced 
in the relative motion of two adjacent fluid layers. 
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Fig. 3. Hydrogen column density (on a logarithmic scale) after 76 000 years for the slow, 
32 km s""*^ [top] and fast, 73 kms""*^ [bottom] models seen at different inclination angles, 
i. Increasing the inclination angle reduces the density contrast between the peak at the 
apex of the bow shock and the rest of the cometary structure, making the detection and 
identification of highly inclined systems more difficult. 



By contrast, the 'fast' model, with nn = 0.3 cm~^ and = 73kms~^, is dom- 
inated by K-H instabilities; the greater stellar motion increases the shear between 
the ISM and stellar wind, reducing the K-H growth timescale. The initially lower 
ISM density results in less cooling and thus less compression; it takes much longer 
to grow R-T instabilities and any excitations that do develop are quickly advected 
downstream. The 'gentle' fluctuations of the K-H instability results in a more 
layered, filamentary appearance (see Fig. [s] [bottom]). 

The appearance of the bow shock is also dependent on the emitting species (see 
Figs. 13 and 14, Mohamed et al. '2012). In our models, the total emissivity is a 
sum of the contributions from 15 different coolants, e.g., rotational and vibrational 
transitions of H2, CO and H2O, H2 dissociative cooling and reformation heating, 
gas-grain cooling/heating, and atomic lines. While some species radiate from 
the entire bow shock surface, e.g., H2O, others are almost entirely confined to 
the reverse shock or forward shock, e.g., CO. Emission primarily from a forward 
shock (hotter gas) results in a much smoother bow shock shell, e.g., the atomic 
line radiation, whereas emission from the reverse shock produces a more layered 
structure, e.g., collisional excitation of H2O with H2. Several coolants, such as 
gas-grain, rotational transitions of CO and H2O, and the heating species produce 
a more 'finger-like', clumpy bow shock sub-structure. 
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Fig. 4. The minimum bow shock shell mass as a function of time for a slow model 
(points) and the fast model (lines). For both models, the RSG and ISM contributions 
are in red and blue, respectively, and their combined mass is plotted in black. The 
masses derived from the IRAS and A i^^Ai?/ observations are ^O.OSM© and ^O.OOSM©, 
respectively (indicated by arrows). 

4 Discussion 

Assuming no other energy sources are present and that during the collision all 
the kinetic energy of both the ISM and stellar wind is thermalised, the theo- 
retical upper limit for the bolometric luminosity of the bow shock is given by 
^tot ^ ^^w'^w + ^^w'^*- For the parameters adopted in our models, £^tot does 
not exceed ~6 x 10^^ ergs s~^. In reality, however, only a small fraction of the 
kinetic energy is radiated from the bow shock; ~16% and ~29% for the fast and 
slow models, respectively. The AKARI (65/im) and the IRAS (60/im) luminosities 
are '^7x10^^ ergss"^ and '^5x10^^ ergss"^, respectively. The latter exceeds the 
theoretical upper limit for the bolometric luminosity by almost an order of magni- 
tude, but is likely overestimated due to contamination from the bar and Betelgeuse 
itself, which is very luminous in the infrared. Whereas the more recent Herschel 
observations are in good agreement with the AKARI values (Decin et al. I2012|) . 
Although the luminosities based on these higher resolution observations are con- 
sistent with the theoretical upper limit, as discussed above, only a small fraction 
of the kinetic energy is thermalised. Furthermore, an even smaller fraction of this 
will be radiated in the far-infrared; from our simulations the combined luminos- 
ity from species thought to be responsible for the far-infrared emission, i.e. dust 
grains, C and O fine structure lines, is at least three orders of magnitude lower 
than the observed flux. The most likely explanation is that Betelgeuse's radiation, 
and the radiation produced by hot gas in the bow shock itself, are absorbed and 
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reemitted by the gas and dust in the far-infrared. 

The evolution of the bow shock shell mass is shown in Fig. [4j Only shocked 
material with x < in the bow shock head is included (recall, in the models the 
star is stationary and positioned at x^y^z= (0,0,0) with the ISM moving in the 
direction of the -\-x axis); this yields a lower limit for the mass in the bow shock 
at any particular time. Assuming the average mass-loss rate of the star has not 
varied significantly, we can compare the bow shock shell mass derived from the 
models with observational estimates and constrain the age of the bow shock. The 
bow shock mass derived from the IRAS fiux is 0.033 Mq and corresponds to a bow 
shock age of ^35 000 years (see Fig. [4j arrows). However, as discussed above, the 
IRAS fiux is likely an overestimate, thus the age is an upper limit. The AKARI 
and Herschel how shock masses are an order of magnitude lower, ~O.OO33M0 
and ~O.OO24M0 (Decin et at. |20T2j), respectively, which would imply an age of 
~10 000 years. If this is the case, however, the wind would not have had sufiicient 
time to expand to the current bow shock radius. One possible solution is that 
the observed shell mass is underestimated due to uncertainties in the fiux to mass 
conversion (e.g., the distance to the star, the gas-to-dust ratio, dust emissivity). 
In our models the shell takes ^20 000 years to reach the observed bow shock 
radius by which time the mass in the bow shock is approximately 0.02 M©. This 
is higher than the value obtained from the far-infrared observations, but may be 
consistent within the uncertainties, and is in agreement with the masses based on 
21cm neutral H observations (Le Bertre et al. 12012) Decin et al. I2012p . Note, 
however, that at ~20 000 years none of our models are close to reaching a steady 
state. 

The shape of Betelgeuse's bow shock is more circular than parabolic. As shown 
m Fig.[3]the bow shock becomes increasingly circular with larger inclination angles, 
i.e. at large angles between the apex of the bow shock and the plane of the sky. 
Ueta et al. (12008^ derived an inclination of 56° using Eqs.[LT]and[L2l(i.e. assumme; 
a steady state) which is consistent with ~50° based on the tangential and radial 
velocities. However, an alternative explanation for the circular shape is that the 
bow shock has not yet reached a steady state. From the Herschel and AKARI 
observations, the ratio of i?(0°)/i?(90°) is approximately 0.7, which is much greater 
than the equilibrium value and corresponds to an age of <30 000 years (Fig. [T]). 

The multiple arcs and even their bright knots in the Herschel observations 
closely resemble the filamentary structure of the fast model (Fig. |2jb]). The fil- 
aments arise when the K-H instabilities are seen projected onto the plane of the 
sky. This projection effect could account for some of the arcs and other struc- 
tures observed at locations well inside the bow shock radius (although the very 
large mass of the HI detached shell would be difficult to explain). From the radial 
velocity and proper motion, the space velocity of Betelgeuse is unlikely to be as 
high as 73 kms~^. (Understanding the origin of the far-ultraviolet emission may 
put constraints on the upper end of the stellar velocity.) The similarity between 
the observations and the fast model, and the lack of clumpy sub-structure that 
characterised the slow models, suggests that the bow shock is dominated by K-H 
rather than R-T instabilities. This situation could occur for the slow models if 
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the RSG wind expands into a much lower density, hot ISM. These conditions re- 
duce the cooling and hence compression in the bow shock, increasing the growth 
time for R-T instabilities. However, for the slow model, such conditions are not 
consistent with the relation = 40 Hj^^^^ kms~^, i.e. the ram pressures are not 
in equilibrium and the bow shock is not yet in a steady state. Thus the overall 
smooth appearance and lack of well-developed instabilities further strengthen the 
argument that the bow shock must be young. 

5 Implications 

We find that many of the physical and morphological characteristics of Betelgeuse's 
bow shock, e.g., the smooth circular shape, the low shell mass and multiple arc sub- 
structure, are consistent with a young bow shock (<30 000 years). Consequently, 
within this time frame, the local ISM through which the star is moving and/or the 
stellar wind must have undergone significant changes. In Mohamed et al. (120121) 
we proposed that such dramatic changes may have occurred if Betelgeuse only 
recently became a RSG, transitioning from either a main sequence (MS) star or a 
blue supergiant (BSC) (i.e. moving from the 'blue' to the 'red' in the Hertzsprung- 
Russell diagram). The radius of Betelgeuse's MS or BSG wind bubble would have 
been ~1 pc, assuming typical wind mass-loss rates (~10~^ M© yr~^) and wind 
velocities (^lO^kms"^) for such hot, blue stars. A RSG phase of a ^fewxlOOOO 
years would bring the star close to the edge of such a bubble; thus, the mysterious 
'bar' ahead of Betelgeuse's bow shock could be a remnant shell produced during 
this earlier phase of 'blue' evolution (Mohamed et al. 2012). A blue-red transition 
would also mean that the RSG wind expands into a MS or BSG bubble filled 
with low density, hot gas, precisely the conditions required to form the observed 
smooth, K-H dominated, multiple arc characteristics of Betelgeuse's bow shock. 

Mackey et al. (j2012p carried out a detailed investigation of a BSG to RSG 
transition, including an evolving wind with a non-constant mass-loss rate and 
wind velocity. Their models reproduce the observed bow shock mass and multiple 
arcs. They also find that the receding BSG bow shock is a plausible candidate 
for the linear bar structure. More detailed comparisons of the models with the 
observations will require a more sophisticated treatment of dust, radiative transfer 
and possibly magnetic fields (see Decin et al. 12012^ . 

Further observations are also required to reduce the number of free and un- 
certain parameters in the models. In particular, a more accurate distance to 
Betelgeuse (197± 45 pc, see Harper et al. I2008P would reduce the uncertainty in 
several key areas, e.g., in deriving the space velocity of the star (to this end a more 
accurate proper motion and radial velocity are also needed). The bow shock mass 
also depends on the distance as well as the highly uncertain dust properties, e.g., 
composition, the dust-to-gas ratio and the dust temperature. Future observations, 
e.g., with ALMA^ could constrain the gas density, temperature and velocity struc- 
ture in the CSE. Indeed, tracing material from the stellar photosphere all the way 
to the bow shock would give us insight into the mass- loss history of Betelgeuse, a 
key ingredient in stellar evolution models. 
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